Package loading

Before performing the analysis, several R-based packages should be installed on your local machine. Note that the packages may require an additional installation of other dependencies. Make sure that you have installed recommended version of the packages and your machine meets the mimimal system requirements for carrying out this analysis.

library(scater)
library(M3Drop)
library(monocle)
library(Seurat)
library(mclust)
library(scran)
library(SC3)

Data information

Zeisel2015 is a UMI-based scRNAseq dataset published in Zeisel et al. (2015) (Link to study: http://science.sciencemag.org/content/347/6226/1138) that contains the expression matrix of endogenous, spike-in and mitochondrial genes from 3005 cells from mouse cortex and hippocampus tissues. Additional metadata information are available, including cells origin, type or tissue. For the purpose of this tutorial we stored Zeisel2015 count matrix with metadata information in two separate txt files. We can read these files and create an easy-to-work SingleCellExperiment experiment object.

all.counts <- read.table(file ="data/all.counts.txt")
metadata <- read.table(file ="data/metadata.txt")

sce <- SingleCellExperiment(list(counts=as.matrix(all.counts)), colData=metadata)
sce
## class: SingleCellExperiment 
## dim: 19896 3005 
## metadata(0):
## assays(1): counts
## rownames(19896): Tspan12 Tshz1 ... ERCC-00170 ERCC-00171
## rowData names(0):
## colnames(3005): V3 V4 ... V3006 V3007
## colData names(10): tissue group.. ... level1class level2class
## reducedDimNames(0):
## spikeNames(0):

We can explore the information about cells and genes stored in SingleCellExperiment experiment object under colData and rawData slots, respectively.

counts(sce)[1:4,1:4]
##          V3 V4 V5 V6
## Tspan12   0  0  0  3
## Tshz1     3  1  0  2
## Fnbp1l    3  1  6  4
## Adamts15  0  0  0  0
names(colData(sce)) #available metadata information
##  [1] "tissue"         "group.."        "total.mRNA.mol" "well"          
##  [5] "sex"            "age"            "diameter"       "cell_id"       
##  [9] "level1class"    "level2class"
table(colData(sce)$tissue) #tissue annotation
## 
## ca1hippocampus       sscortex 
##           1314           1691
table(colData(sce)$level1class) #cell type annotation
## 
## astrocytes_ependymal    endothelial-mural         interneurons 
##                  224                  235                  290 
##            microglia     oligodendrocytes        pyramidal CA1 
##                   98                  820                  939 
##         pyramidal SS 
##                  399
names(rowData(sce))
## character(0)

As mentioned, Zeisel2015 dataset contains also the expression of spike-ins and mitochondrial genes. Spike-in gene names usually starts from “ERCC” and mitochondrial genes contain “mt” string.

is.spike <- grepl("^ERCC", rownames(sce))
summary(is.spike)
##    Mode   FALSE    TRUE 
## logical   19839      57
is.mito <- grepl("^mt-", rownames(sce))
summary(is.mito)
##    Mode   FALSE    TRUE 
## logical   19862      34

Quality control and data preprocessing

We can use scater package to automatically calculate several quality metrics per cell and per gene. These metrics will be useful to i.e. detect low quality cells or lowly expressed genes that should be removed from further analysis.

sce <- calculateQCMetrics(sce,  feature_controls=list(Spike=is.spike, Mt=is.mito))
colnames(colData(sce))
##  [1] "tissue"                                        
##  [2] "group.."                                       
##  [3] "total.mRNA.mol"                                
##  [4] "well"                                          
##  [5] "sex"                                           
##  [6] "age"                                           
##  [7] "diameter"                                      
##  [8] "cell_id"                                       
##  [9] "level1class"                                   
## [10] "level2class"                                   
## [11] "is_cell_control"                               
## [12] "total_features_by_counts"                      
## [13] "log10_total_features_by_counts"                
## [14] "total_counts"                                  
## [15] "log10_total_counts"                            
## [16] "pct_counts_in_top_50_features"                 
## [17] "pct_counts_in_top_100_features"                
## [18] "pct_counts_in_top_200_features"                
## [19] "pct_counts_in_top_500_features"                
## [20] "total_features_by_counts_endogenous"           
## [21] "log10_total_features_by_counts_endogenous"     
## [22] "total_counts_endogenous"                       
## [23] "log10_total_counts_endogenous"                 
## [24] "pct_counts_endogenous"                         
## [25] "pct_counts_in_top_50_features_endogenous"      
## [26] "pct_counts_in_top_100_features_endogenous"     
## [27] "pct_counts_in_top_200_features_endogenous"     
## [28] "pct_counts_in_top_500_features_endogenous"     
## [29] "total_features_by_counts_feature_control"      
## [30] "log10_total_features_by_counts_feature_control"
## [31] "total_counts_feature_control"                  
## [32] "log10_total_counts_feature_control"            
## [33] "pct_counts_feature_control"                    
## [34] "pct_counts_in_top_50_features_feature_control" 
## [35] "pct_counts_in_top_100_features_feature_control"
## [36] "pct_counts_in_top_200_features_feature_control"
## [37] "pct_counts_in_top_500_features_feature_control"
## [38] "total_features_by_counts_Spike"                
## [39] "log10_total_features_by_counts_Spike"          
## [40] "total_counts_Spike"                            
## [41] "log10_total_counts_Spike"                      
## [42] "pct_counts_Spike"                              
## [43] "pct_counts_in_top_50_features_Spike"           
## [44] "pct_counts_in_top_100_features_Spike"          
## [45] "pct_counts_in_top_200_features_Spike"          
## [46] "pct_counts_in_top_500_features_Spike"          
## [47] "total_features_by_counts_Mt"                   
## [48] "log10_total_features_by_counts_Mt"             
## [49] "total_counts_Mt"                               
## [50] "log10_total_counts_Mt"                         
## [51] "pct_counts_Mt"                                 
## [52] "pct_counts_in_top_50_features_Mt"              
## [53] "pct_counts_in_top_100_features_Mt"             
## [54] "pct_counts_in_top_200_features_Mt"             
## [55] "pct_counts_in_top_500_features_Mt"

Quality metrics per cell

colData(sce)$total_counts[1:5] #Total count for each cell
## [1] 29060 29304 38929 40557 27625
colData(sce)$total_features_by_counts[1:5] #Number of genes with non zero counts per cell
## [1] 4898 4750 6090 5887 4753
colData(sce)$pct_counts_Spike[1:5] #Percentage of spike-in counts per cell
## [1] 23.07639 21.95946 16.27322 17.33856 21.46968

Quality metrics per gene

rowData(sce)$mean_counts[1:5] #Average count per gene
## [1] 0.27886855 0.42362729 1.04292845 0.04459235 0.37936772
rowData(sce)$n_cells_by_counts[1:5] #Number of cells with non zero counts per gene
## [1]  472  573 1167   77  669
rowData(sce)$pct_dropout_by_counts[1:5] #Percentage of dropouts per gene
## [1] 84.29285 80.93178 61.16473 97.43760 77.73710

Plot histograms based on quality metrics

One can plot histograms based on quality metrics to determine possible thresholds for cell/gene filterings.

hist(sce$total_counts, breaks=50,  xlab="Library size", ylab="Number of cells")

hist(sce$total_features_by_counts, xlab="Number of expressed genes", breaks=50,  ylab="Number of cells")

hist(sce$pct_counts_Spike, xlab="ERCC proportion (%)", ylab="Number of cells", breaks=50)

Automaticly filter outlier cells based on quality metrics

Alternatively, isOutlier function provides an automatic detection of low quality cells based on quality metrics. In this example, we are filtering cells based on library size, number of expressed genes per cell and total count over all spike-in transcripts in each cell. Excluded cells are those with the total number of expressed genes and the total sum of counts more than 3 median absolute deviations below the median across the genes.

libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
spike.drop <- isOutlier(sce$pct_counts_Spike, nmads=3, type="higher")
sce <- sce[,!(libsize.drop | feature.drop | spike.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop), BySpike=sum(spike.drop), Remaining=ncol(sce))
##   ByLibSize ByFeature BySpike Remaining
## 1         8         3       8      2989

Filter spike-in genes

After quality control on cells, we can remove the expression of spike-in genes from the count matrix.

sce=sce[-which(grepl("^ERCC-", rownames(sce))),]
sce
## class: SingleCellExperiment 
## dim: 19839 2989 
## metadata(0):
## assays(1): counts
## rownames(19839): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
## rowData names(9): is_feature_control is_feature_control_Spike ...
##   total_counts log10_total_counts
## colnames(2989): V3 V4 ... V3006 V3007
## colData names(55): tissue group.. ...
##   pct_counts_in_top_200_features_Mt
##   pct_counts_in_top_500_features_Mt
## reducedDimNames(0):
## spikeNames(0):

Filter lowly expressed genes based on the average count per gene

Now, we should filter out lowly expressed genes as they do not provide any insights into the underlying biology. In the filtering we removed lowly expressed genes that are genes with average expression count (adjusted by library size) equal to 0.

rowData(sce)$ave.counts <- calcAverage(sce, exprs_values = "counts", use_size_factors=FALSE)
to.keep <- rowData(sce)$ave.counts > 0
sce <- sce[to.keep,]
summary(to.keep)
##    Mode   FALSE    TRUE 
## logical       2   19837
dim(sce)
## [1] 19837  2989

Explore highly expressed genes, if required filter overexpressed features

We can also explore highly and overexpressed genes. The overexpressed genes, are those genes which expression is higher than others by several magnitudes. Overexpressed genes can bias further procedures such as clustering, therefore, they should be removed before performing downstream analysis. Note that one of the most common overexpressed gene is “Malat1”.

plotHighestExprs(sce, exprs_values = "counts")

sce <- sce[-which(rownames(sce)=="Malat1"),]
dim(sce)
## [1] 19836  2989

Perform standard normalization by computing counts per million (CPM)

After quality control and filtering, one should normalize the dataset to remove potential technical bias. One strategy is to use CPM (Count Per Million) normalization which is a correction to remove the noise related to sequencing depth. CPM divides each count by its total sum (across all the genes) and multiplies by one million. In this way, each cell has the same total sum of the counts.

cpm(sce) <- calculateCPM(sce, use_size_factors=FALSE)
cpm(sce)[1:4,1:4]
##                V3       V4       V5        V6
## Tspan12    0.0000  0.00000   0.0000  92.95696
## Tshz1    139.3275 45.45455   0.0000  61.97131
## Fnbp1l   139.3275 45.45455 191.3448 123.94261
## Adamts15   0.0000  0.00000   0.0000   0.00000

Perform normalization using decnvolution method

Another strategy, more suitable for scRNAseq data, is deconvolution method implemented in scran package. The deconvolution method normalizes data by cells-pooled size factors that account for dropout biases. More details about this normalization technique can be found in study https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0947-7.

set.seed(100)
clusters <- quickCluster(sce, min.size=200, min.mean=0.1, method="igraph")
sce <- computeSumFactors(sce, cluster=clusters, min.mean=0.1)
sce <- scater::normalize(sce, return_log = FALSE)
normcounts(sce)[1:4,1:4]
##                V3        V4       V5        V6
## Tspan12  0.000000 0.0000000 0.000000 1.2865579
## Tshz1    1.896026 0.6416602 0.000000 0.8577053
## Fnbp1l   1.896026 0.6416602 2.532498 1.7154106
## Adamts15 0.000000 0.0000000 0.000000 0.0000000

Perform feature selection to retrieve most variable genes

Feature selection step is aimed at preserving biologically relevant information in the dataset and improving computational efficiency of downstream analyses. In most of the cases, it seeks for highly variable genes based on mean/variance relationship. M3Drop is one of the packages for automatic detection of most variable genes. It first calculates the mean and square coefficient of variation and fits the quadratic curve between the two variables. Then a chi-square test is used to find high variable genes which are significantly above the curve.

hvg_genes <- BrenneckeGetVariableGenes(normcounts(sce), fdr = 0.01, minBiolDisp = 0.5)

length(hvg_genes)
## [1] 14225

scran package also allows to identify genes that drive biological heterogeneity by modelling per-gene variance. The method first decomposes the total variance of each gene into its biological and technical components then fits mean-variance trend to the normalized log-expression values. By default, the method uses spike-in genes to identify technical noise. If no spike-ins are available, the method will use all the features with assumption that technical noise is the major contributor to the variance of most genes in the dataset. Once technical component is estimated it will be substracted from the total variance to obtain the biological component for each gene. Genes with high biological components are of interest. Note that the method requires to input log-transformed normalized counts and outputs an ordered list of all genes from which the user has to extract a certain number of top features for use in downstream analysis.

assay(sce, "logcounts") <- log2(normcounts(sce)+1) 

fit <- trendVar(sce, parametric=TRUE, use.spikes=FALSE)
decomp <- decomposeVar(sce, fit)
head(decomp)
## DataFrame with 6 rows and 6 columns
##                        mean              total                  bio
##                   <numeric>          <numeric>            <numeric>
## Tspan12    0.21954217578161  0.367997453453532   0.0676731873287152
## Tshz1     0.311930968051195  0.549001316751025    0.138350625999446
## Fnbp1l     0.61796543640819  0.839797618303715   0.0845728838805125
## Adamts15  0.035777983570672 0.0699552575563804   0.0165716888849863
## Cldn12    0.279720893210915  0.383207030299019    0.010109269426282
## Rxfp1    0.0365498687107355 0.0473109676070912 -0.00722212286053495
##                        tech              p.value                  FDR
##                   <numeric>            <numeric>            <numeric>
## Tspan12   0.300324266124817 2.00236966401108e-16 1.16889360374702e-15
## Tshz1     0.410650690751579 1.90913067042286e-32 1.58915299951775e-31
## Fnbp1l    0.755224734423203 1.42457379368776e-05 5.61786198242353e-05
## Adamts15 0.0533835686713941 3.28461490021999e-28  2.5107368462722e-27
## Cldn12    0.373097760872737    0.147643066563371    0.420420308405258
## Rxfp1    0.0545330904676262    0.999999957325667                    1
genes_ordered <- order(decomp$bio, decreasing=TRUE)
hvg_genes <- genes_ordered[1:500]
length(hvg_genes)
## [1] 500

Data visualization

For reducing dimension of the data you can use PCA, tSNE or UMAP techniques implemented in scater package. Note that projections should be applied to normalized and log-transformed count matrices. For reproducibility of tSNE and UMAP projections - set the seed for generating random variables.

plotPCA(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))  
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

set.seed(100)
plotTSNE(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

plotUMAP(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

Exercise

How the projections would change if you would plot raw or CPM counts followed by log-transformation?

Additional plots

You can use colour_by argument to color points in any of the above projections by a selected marker gene or any information stored in the annotation. Additional violin plots allow you to i.e. visualize the expression of the selected gene across all the cells from a given cell type. For the illustrative purpose, we used Gad1 gene taken from the literature, which is a marker for interneuron cell type in the cortex tissue.

set.seed(100)
plotTSNE(sce, colour_by="Gad1", run_args=c(exprs_values = "logcounts")) #marker for interneurons
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

plotExpression(sce, "Gad1", x = "level1class", colour_by = "level1class", exprs_values = "logcounts") 

Clustering

To cluster cells we can use SC3 package. SC3 is a consensus clustering method for single-cell RNA-seq data. It first calculates distances between the cells using three metrics: Euclidean, Pearson, and Spearman. On each of the obtained distance matrix two transformations are applied: PCA and Laplacian graph. Then K-Means clustering technique is used to cluster the transformed distance matrices subject to the first d eigenvectors. In result, several individual clusterings are obtained which are further combined into a single consensus clustering using Cluster-based Similarity Partitioning Algorithm (CSPA). Note that when using SC3, one can first estimate the number of populations to be used in the clustering.

sce <- sc3_estimate_k(sce)
## Estimating k...
k_input <- metadata(sce)$sc3$k_estimation
k_input #number of estimated clusters
## [1] 32
rowData(sce)$feature_symbol <- rownames(sce)
sce <- sc3(sce, ks = k_input, gene_filter = FALSE, biology = TRUE)
## Setting SC3 parameters...
## Your dataset contains more than 2000 cells. Adjusting the nstart parameter of kmeans to 50 for faster performance...
## Calculating distances between the cells...
## Performing transformations and calculating eigenvectors...
## Performing k-means clustering...
## Calculating consensus matrix...
## Calculating biology...
table(colData(sce)[paste0("sc3_",k_input, "_clusters")][[1]])
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
## 335  66 286  11   2  26  13  92  86  30   6  80  92 103   6  10   9 182 
##  19  20  21  22  23  24  25  26  27  28  29  30  31  32 
##   9  71  80 226  41  64  56 135   4 111 304  32  43 378

Evaluate clustering output based on annotation

If the cell annotation is available, one can measure the effectiveness of clustering output using Adjusted Rand Index (ARI). ARI index measures similarity between the partition obtained from clustering and the partition obtained from dataset annotation. The values of the ARI range can be negative if the agreement of the partitions is worse then the agreement expected by chance, or between 0 and 1 for clustering better then chance. ARI close to 1 indicate high accuracy of the method in detecting annotated cell populations.

adjustedRandIndex(colData(sce)[paste0("sc3_",k_input, "_clusters")][[1]],sce$level1class)
## [1] 0.3580484

Exercise

How the clustering would change if you would use true number of cell populations (taken from annotation) in the inference of SC3?

Complete data analysis using Monocle package

Create CellDataSet data object

CellDataSet data object is similar for SingleCellExperiment object. It stores information about the experiment (i.e. expression matrix and metadata) in “slots”. Note that when creating a new CellDataSet data should not be normalized - Monocle, when performing downstream analysis steps, normalizes the data internally.

rowData(sce)$gene_short_name <- rownames(sce)
pd <- new("AnnotatedDataFrame", data = data.frame(colData(sce))) #cell info
fd <- new("AnnotatedDataFrame", data = data.frame(rowData(sce))) #gene info
cds <- newCellDataSet(cellData = counts(sce), phenoData = pd, featureData = fd, expressionFamily=VGAM::negbinomial.size())  
cds
## CellDataSet (storageMode: environment)
## assayData: 19836 features, 2989 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: V3 V4 ... V3007 (2989 total)
##   varLabels: tissue group.. ... Size_Factor (58 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: Tspan12 Tshz1 ... mt-Nd4l (19836 total)
##   fvarLabels: is_feature_control is_feature_control_Spike ...
##     gene_short_name (17 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

Extract highly variable genes for further analysis

cds <- BiocGenerics::estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
## Removing 132 outliers
disp <- dispersionTable(cds)
head(disp)
##    gene_id mean_expression dispersion_fit dispersion_empirical
## 1  Tspan12      0.27371308       5.352092             9.664429
## 2    Tshz1      0.41233507       3.788779             8.705790
## 3   Fnbp1l      0.83029734       2.234915             2.957010
## 4 Adamts15      0.04653374      28.054145            71.442708
## 5   Cldn12      0.30721173       4.845039             5.144358
## 6    Rxfp1      0.03723447      34.885324            14.381728
hvg <- subset(disp, mean_expression >= 0.5 & dispersion_empirical >= 0.1)
cds <- setOrderingFilter(cds, ordering_genes = hvg$gene_id)
plot_ordering_genes(cds)
## Warning: Transformation introduced infinite values in continuous y-axis

cds <- cds[hvg$gene_id,]

Reduce dimension of the data and cluster cells to detect cell populations

cds <- reduceDimension(cds, max_components = 2, reduction_method = 'tSNE', verbose = T)
## Warning in if (method == "PCA") {: the condition has length > 1 and only
## the first element will be used
cds <- clusterCells(cds, num_clusters = NULL)
## Distance cutoff calculated to 3.325859

Compare clusterization with annotated cell types

p1 <- monocle::plot_cell_clusters(cds, color_by = 'Cluster')
p1

p2 <- monocle::plot_cell_clusters(cds, color_by = 'level1class')
p2

Compute Adjusted Rand Index

adjustedRandIndex(cds$Cluster,cds$level1class)
## [1] 0.7575202

Complete data analysis using Seurat package

Create Seurat object

seur <- CreateSeuratObject(raw.data = counts(sce), meta.data = data.frame(colData(sce)), min.cells = 0, min.genes = 0)
seur
## An object of class seurat in project SeuratProject 
##  19836 genes across 2989 samples.

Seurat require to log-normalize and scale data before performing clustering

seur <- NormalizeData(seur, normalization.method = "LogNormalize")
seur@data[1:4,1:4]
## 4 x 4 sparse Matrix of class "dgCMatrix"
##                 V3        V4       V5        V6
## Tspan12  .         .         .        0.6572970
## Tshz1    0.8726628 0.3746934 .        0.4822490
## Fnbp1l   0.8726628 0.3746934 1.069337 0.8062196
## Adamts15 .         .         .        .

Extract highly variable genes for further analysis

seur <- FindVariableGenes(object = seur)

Scale the data

seur<- Seurat::ScaleData(object = seur)
## Scaling data matrix

Reduce dataset dimension with PCA and clusterize cells

seur <- RunPCA(object = seur, pc.genes=seur@var.genes)
seur <- FindClusters(object = seur, reduction.type="pca")
table(seur@ident)
## 
##   0   1   2   3   4   5   6 
## 848 841 366 318 276 257  83

Compare clusterization with annotated cell types

Seurat::DimPlot(seur, group.by = "ident", reduction.use = "pca")

Seurat::DimPlot(seur, group.by = "level1class", reduction.use = "pca")

Compute Adjusted Rand Index

adjustedRandIndex(as.numeric(seur@ident), seur@meta.data$level1class)
## [1] 0.8072137

Package versions used in this analysis

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
## 
## Matrix products: default
## BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] SC3_1.10.1                  scran_1.10.2               
##  [3] mclust_5.4.3                Seurat_2.3.4               
##  [5] cowplot_0.9.4               monocle_2.99.2             
##  [7] L1Graph_0.1.1               lpSolveAPI_5.5.2.0-17      
##  [9] DDRTree_0.1.5               irlba_2.3.3                
## [11] igraph_1.2.4                Matrix_1.2-17              
## [13] M3Drop_1.8.1                numDeriv_2016.8-1          
## [15] scater_1.10.1               ggplot2_3.1.0              
## [17] SingleCellExperiment_1.4.1  SummarizedExperiment_1.12.0
## [19] DelayedArray_0.8.0          BiocParallel_1.16.6        
## [21] matrixStats_0.54.0          Biobase_2.42.0             
## [23] GenomicRanges_1.34.0        GenomeInfoDb_1.18.2        
## [25] IRanges_2.16.0              S4Vectors_0.20.1           
## [27] BiocGenerics_0.28.0        
## 
## loaded via a namespace (and not attached):
##   [1] prabclus_2.2-7           R.methodsS3_1.7.1       
##   [3] coda_0.19-2              pkgmaker_0.27           
##   [5] tidyr_0.8.3              acepack_1.4.1           
##   [7] bit64_0.9-7              knitr_1.22              
##   [9] R.utils_2.8.0            data.table_1.12.0       
##  [11] rpart_4.1-13             RCurl_1.95-4.12         
##  [13] doParallel_1.0.14        metap_1.1               
##  [15] snow_0.4-3               RANN_2.6.1              
##  [17] VGAM_1.1-1               combinat_0.0-8          
##  [19] proxy_0.4-23             future_1.12.0           
##  [21] bit_1.1-14               webshot_0.5.1           
##  [23] httpuv_1.5.0             assertthat_0.2.1        
##  [25] viridis_0.5.1            xfun_0.5                
##  [27] evaluate_0.13            promises_1.0.1          
##  [29] DEoptimR_1.0-8           caTools_1.17.1.2        
##  [31] DBI_1.0.0                htmlwidgets_1.3         
##  [33] sparsesvd_0.1-4          spdep_1.0-2             
##  [35] purrr_0.3.2              RSpectra_0.13-1         
##  [37] crosstalk_1.0.0          dplyr_0.8.0.1           
##  [39] backports_1.1.3          trimcluster_0.1-2.1     
##  [41] gbRd_0.4-11              deldir_0.1-16           
##  [43] Cairo_1.5-10             ROCR_1.0-7              
##  [45] withr_2.1.2              robustbase_0.93-4       
##  [47] checkmate_1.9.1          cluster_2.0.7-1         
##  [49] ape_5.3                  segmented_0.5-3.0       
##  [51] lazyeval_0.2.2           crayon_1.3.4            
##  [53] hdf5r_1.1.1              labeling_0.3            
##  [55] glmnet_2.0-16            edgeR_3.24.3            
##  [57] pkgconfig_2.0.2          slam_0.1-45             
##  [59] units_0.6-2              nlme_3.1-137            
##  [61] vipor_0.4.5              nnet_7.3-12             
##  [63] rlang_0.4.0              globals_0.12.4          
##  [65] diptest_0.75-7           miniUI_0.1.1.1          
##  [67] registry_0.5-1           doSNOW_1.0.16           
##  [69] ggrastr_0.1.7            lmtest_0.9-36           
##  [71] rngtools_1.3.1           Rhdf5lib_1.4.3          
##  [73] boot_1.3-20              zoo_1.8-5               
##  [75] base64enc_0.1-3          beeswarm_0.2.3          
##  [77] ggridges_0.5.1           pheatmap_1.0.12         
##  [79] png_0.1-7                viridisLite_0.3.0       
##  [81] bitops_1.0-6             R.oo_1.22.0             
##  [83] KernSmooth_2.23-15       DelayedMatrixStats_1.4.0
##  [85] rgl_0.100.19             doRNG_1.7.1             
##  [87] classInt_0.3-1           lars_1.2                
##  [89] stringr_1.4.0            manipulateWidget_0.10.0 
##  [91] scales_1.0.0             magrittr_1.5            
##  [93] plyr_1.8.4               ica_1.0-2               
##  [95] gplots_3.0.1.1           bibtex_0.4.2            
##  [97] gdata_2.18.0             zlibbioc_1.28.0         
##  [99] compiler_3.5.1           HSMMSingleCell_1.2.0    
## [101] lsei_1.2-0               bbmle_1.0.20            
## [103] RColorBrewer_1.1-2       rrcov_1.4-7             
## [105] fitdistrplus_1.0-14      dtw_1.20-1              
## [107] XVector_0.22.0           LearnBayes_2.15.1       
## [109] listenv_0.7.0            pbapply_1.4-0           
## [111] htmlTable_1.13.1         Formula_1.2-3           
## [113] MASS_7.3-51.3            tidyselect_0.2.5        
## [115] stringi_1.4.3            densityClust_0.3        
## [117] yaml_2.2.0               locfit_1.5-9.1          
## [119] latticeExtra_0.6-28      ggrepel_0.8.0           
## [121] pbmcapply_1.3.1          grid_3.5.1              
## [123] tools_3.5.1              rstudioapi_0.10         
## [125] foreach_1.4.4            foreign_0.8-71          
## [127] gridExtra_2.3            Rtsne_0.15              
## [129] digest_0.6.18            FNN_1.1.3               
## [131] shiny_1.2.0              qlcMatrix_0.9.7         
## [133] fpc_2.1-11.1             Rcpp_1.0.1              
## [135] SDMTools_1.1-221         later_0.8.0             
## [137] WriteXLS_4.1.0           httr_1.4.0              
## [139] sf_0.7-3                 npsurv_0.4-0            
## [141] kernlab_0.9-27           Rdpack_0.10-1           
## [143] colorspace_1.4-1         reticulate_1.11.1       
## [145] umap_0.2.0.0             splines_3.5.1           
## [147] statmod_1.4.30           expm_0.999-4            
## [149] sp_1.3-1                 flexmix_2.3-15          
## [151] plotly_4.8.0             spData_0.3.0            
## [153] xtable_1.8-3             jsonlite_1.6            
## [155] dynamicTreeCut_1.63-1    modeltools_0.2-22       
## [157] R6_2.4.0                 gmodels_2.18.1          
## [159] Hmisc_4.2-0              pillar_1.4.2            
## [161] htmltools_0.3.6          mime_0.6                
## [163] glue_1.3.1               BiocNeighbors_1.0.0     
## [165] class_7.3-15             codetools_0.2-16        
## [167] pcaPP_1.9-73             tsne_0.1-3              
## [169] mvtnorm_1.0-10           lattice_0.20-38         
## [171] tibble_2.1.1             mixtools_1.1.0          
## [173] ggbeeswarm_0.6.0         gtools_3.8.1            
## [175] survival_2.44-1.1        limma_3.38.3            
## [177] rmarkdown_1.12           docopt_0.6.1            
## [179] fastICA_1.2-1            munsell_0.5.0           
## [181] e1071_1.7-1              rhdf5_2.26.2            
## [183] GenomeInfoDbData_1.2.0   iterators_1.0.10        
## [185] HDF5Array_1.10.1         reshape2_1.4.3          
## [187] gtable_0.3.0